Efficient time series detection of the strong stochasticity threshold in 
Fermi-Pasta-Ulam oscillator lattices 
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In this work we study the possibility of detecting the so-called strong stochasticity threshold, 
i.e. the transition between weak and strong chaos as the energy density of the system is increased, 
in anharmonic oscillator chains by means of the 0-1 test for chaos. We compare the result of the 
aforementioned methodology with the scaling behavior of the largest Lyapunov exponent computed 
by means of tangent space dynamics, that has so far been the most reliable method available to detect 
the strong stochasticity threshold. We find that indeed the 0-1 test can perform the detection in 
the range of energy density values studied. Furthermore, we determined that conventional nonlinear 
time series analysis methods fail to properly compute the largest Lyapounov exponent even for very 
large data sets, whereas the computational effort of the 0-1 test remains the same in the whole 
range of values of the energy density considered with moderate size time series. Therefore, our 
results show that, for a qualitative probing of phase space, the 0-1 test can be an effective tool if its 
limitations are properly taken into account. 

PACS numbers: 05.45.Tp, 05.45.Pq, 05.45.Jn 



I. INTRODUCTION 

Deterministic chaos, defined as a dynamical regime 
with sensitive dependence on initial conditions, can be 
readily quantified by the largest Lyapunov exponent 
(LLE) A, which is an averaged exponential growth rate 
of linear perturbations to the considered motion jlj . The 
LLE is easily computed in numerical simulation, but it is 
difficult to obtain in experiment, where usually there is 
no knowledge of the nonlinear equations that govern the 
time evolution of the system. Hence, a number of alterna- 
tive chaos-detection techniques that rely on the temporal 
evolution of a reduced set of variables have been devel- 
oped over the years with varying degrees of success @j. 

The 0-1 test Q, which uses as input the time series of 
an observed variable, was designed to distinguish regu- 
lar behavior, if is the output, from chaotic dynamics, 
if the output is 1, in deterministic dynamical systems. 
It has been successfully applied to both numerically ob- 
tained as well as to simple experimental data 
However, in a previous work jlpj we have presented ev- 
idence that, when applied to time series stemming from 
nonhomogeneous, i.e dissimilar mass, many-degrees-of- 
freedom systems, the test has serious limitations as an 
efficient chaos detector for a large system size N. On the 
other hand, for small homogeneous systems, but in a very 
low energy regime (indistinguishable of regular dynam- 
ics for very large time scales and termed weak chaos), 
the test misclassifies the provided signal as regular, even 
though the LLE has a vanishingly small non-negative 
value. 
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Of the foregoing results, the latter is by far the most 
important to be taken into account when applying the 0-1 
test, since the homogeneity condition can be readily con- 
trolled in a simulation or in an experiment, and systems 
with a small number of homogeneous components are in- 
deed physically relevant, as will be discussed below. The 
analysis of Ref. [lOj makes it plausible to infer that, for 
homogeneous mass systems composed of a small number 
of particles, it can be generically expected that, in a low 
energy (or temperature) regime, the 0-1 test will render a 
result close to zero notwithstanding the actual dynamical 
regime of the considered system. Now, if a microscopic 
model is available, in principle the LLE provides a way to 
know, with no uncertainty, the type of dynamical regime 
wherein the system is, since the dynamical equations give 
information of all the existing range of scales for arbi- 
trarily long time intervals. However, there are evident 
practical problems with the computation of such quan- 
tity, since the aforementioned conditions, besides being 
unattainable when dealing with experimental data, may 
also result of no physical interest whatsoever, since it is 
now clear that the LLE is not completely satisfactory for 
a proper characterization of the many faces of complexity 
and predictability of many-degrees-of-freedom systems: 
see |11| and references therein, where some examples are 
discussed showing that systems with different dynamics 
can give similar results when analyzed from a time series 
point of view. From this perspective, what can be validly 
ascertained is that, if the result of the 0-1 test is for 
a signal obtained from a simulation or experiment of a 
homogeneous system in a very low energy or tempera- 
ture regime, the dynamics can be regarded as regular for 
the considered time series length, regardless of what the 
outcome of the test would be with a larger (and maybe 
inaccessible) data set. 
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The precise knowledge of the actual dynamical regime 
may not be a crucial piece of information if the main in- 
terest is to probe other phase-space structural character- 
istics of many-degrees-of- freedom systems that go beyond 
the simple distinction between order and chaos. In par- 
ticular, the so called Fermi-Pasta-Ulam (FPU) oscillator 
chain presents a transition from weak to strong chaos as 
the energy density e = E/N of the lattice increases be- 
yond a threshold value known as the strong stochas- 
ticity threshold (SST) [f|[l3j]. The LLE A(e) exhibits a 
change in its scaling behavior precisely at the value e T 
of the SST. If a complete knowledge of the dynamics is 
available, even an analytical estimate of the LLE energy 
dependence for both low and high energies can be read- 
ily obtained fl4j |. However, with only limited information 
(such as the time series of the position and/or momen- 
tum of a single oscillator in an actual experiment) the 
computation of the LLE becomes problematic. There- 
fore, it would be interesting to explore the possibility 
that the 0-1 test could detect the aforementioned transi- 
tion, regardless the misclassification (which, in principle, 
can only be entirely avoided in the infinite time limit) 
of the actual dynamical regime for very low e values ob- 
tained by means of this method. This result could be 
specially relevant for the case of experimental data with 
no a priori information about the microscopic dynamics, 
but with precise structural information available that can 
rule out the limitations addressed in Ref. [Toj . 

This paper is organized as follows. In Sec.[H]we present 
the relevant detail of the 0-1 test and in Sec. IIIII we de- 
scribe the employed model and the relevant details of its 
numerical integration. Sec. IIVI presents the results of the 
detection of the SST by means of the 0-1 test as well as 
a comparison with results stemming from standard non- 
linear time series analysis methods. In Sec. [V] we discuss 
the previous results and present our conclusions. 



II. THE 0-1 TEST FOR CHAOS 

We follow the implementation of the 0-1 test reported 
in Q: See Ref. [l5| for the relevant mathematical re- 
sults that justify the version herein employed. The 
test takes as input a finite data set {(j>(ta)}a=i sam- 
pled at discrete times t a = ar, with sampling time r. 
For a given c s 5ft we construct the transformed series 
£(*a) = S"=i ^(tj) cos(jc), a = 1,2,3,... Now, from the 
mean square displacement defined as 



M(t a )= Km 

7v->oo yv - 
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(1) 



we compute numerically the asymptotic growth rate K = 
lim Q _ i . 00 (logM(t ce ))/logi a by performing a least square 
fit of log M(t a ) versus \ogt a in the range 1 < a < Afi for 
a choice of M such that 1 < M < M and Ni = W/10. 
We compute K for 100 random values of c and the final K 
value is taken as the median of the computed set. Then 



K « stands for regular dynamics and K w 1 for chaotic 
dynamics. 



III. THE MODEL AND ITS NUMERICAL 
SIMULATION 

The Hamiltonian of the one-dimensional N coupled 
FPU lattice from which the employed time series were 
obtained reads as 



N r 



Pi 

2rrij 



(2) 



where {m i ,q i ,p i }fL 1 are the dimensionless mass, dis- 
placement, and momentum of the i oscillator, respec- 
tively. Periodic boundary conditions are assumed. Equi- 
librium displacements {qi{0) = 0} and momenta {pi(0)} 
distributed according to a Maxwell-Boltzmann distribu- 
tion consistent with the desired energy density e are taken 
as initial conditions. The homogeneity of the system is 
assured by taking a unit mass value m, = 1 for all os- 
cillators. Next, the 2N first-order Hamilton equations 
of motion were integrated using a symmetrical version 
of the velocity Verlet algorithm (VVA) [16]. With the 
adopted value At = 0.01 of the time step a faithful rep- 
resentation of a Hamiltonian flow and a driftless value of 
total energy E for the studied time scales is ensured. 



IV. RESULTS 

A. Comparison with tangent space methods 

For each considered e value the position and momen- 
tum time series of a single oscillator of a FPU lattice 
with N = 32 were recorded with a sampling time of 
t = 1 (that corresponds to 100 time steps for At = 0.01) 
with a total time series length of J\f = 10 5 . The 0-1 test 
was applied to the so obtained data sets and the corre- 
sponding asymptotic growth rates K q (position) and K p 
(momentum) for each e was computed. The results are 
displayed in Fig. Ufa). As can be readily appreciated, 
K grows steadily and gradually reaches a ~ 1 value as 
the energy density changes from low to high values. The 
transition between the aforementioned behaviors occurs 
at e G [0.1,0.2]. A further refinement is obtained if one 
observes that K q and K p have different (but albeit very 
close) values up to the threshold e 01 ~ 0.2, where the 
difference K q — K p becomes minimal and K q ,K p « 0.9 
consistently hereafter, as seen in the inset of the same fig- 
ure. The behavior of K vs e indicates a transition from a 
very ordered dynamical regime to a strongly chaotic one, 
with a very precise threshold value e 01 separating both 
regimes. This capability of the 0-1 test to probe detailed 
dynamical information has not been previously acknowl- 
edged. The exact quantification of the degree of expo- 
nential divergence, which also conveys information of the 
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FIG. 1: (Color online) (a) Asymptotic growth rate K vs en- 
ergy density e computed from the position (circles) and mo- 
mentum (squares) time series of a single oscillator of a FPU 
chain with N = 32. The inset presents in more detail the 
transition region where the K values computed from both 
time series become almost indistinguishable. The horizontal 
dashed line indicates the value 0.9 to which both K q and K p 
converge for large e values, (b) LLE vs e obtained by means 
of the standard method. Diamonds correspond to the LLE 
computed using the VVA of Ref . [l6l ] , whereas circles refer to 
data obtained by means of the BSA3 of Ref. [2ll |. Continuous 
lines in (b) refer to the scaling laws valid for low and high e 
values, whereas vertical dot-dashed lines in both (a) and (b) 
indicate the approximate value e T corresponding to the SST. 

aforementioned transition, is obtained by means of the 
LLE. These results for LLE A(e), computed by means of 
the standard method [13, UK] ; are presented in Fig. QJb) . 
The SST e T , obtained from the intersection of the scal- 
ing laws (whose theoretical explanation is considered in 
Ref. [14|) valid for low and high e values, has a remarkable 
agreement with e 01 obtained by means of the 0-1 test and 
is consistent with other estimates [l9j . Since it is known 
that a precise computation of the LLE is difficult in the 
low energy density regime [20l |. we have performed some 
further computations with the more sophisticated third 
order bilateral symplectic algorithm (BSA3) [21|. An ex- 
cellent agreement with the results obtained by means of 
the VVA can be readily appreciated in that same figure. 

B. Results of phase-space reconstruction methods 

Notice that the knowledge of the tangent space dynam- 
ics is not usually available when studying experimental 



time series, since it can only be obtained from an explicit 
model. Therefore, the LLE has to be obtained from dif- 
ferent approximation techniques such as phase space re- 
construction methods, successfully applied to dissipativc 
systems Q. In order to corroborate if the LLE of the 
Hamiltonian model herein studied can indeed be com- 
puted by means of these techniques we proceed to deter- 
mine an embedding dimension to that can approximate 
the underlying dimensionality 2N of the considered dy- 
namical system (22j . In Fig. [5] we present the results 
for both the correlation integral C m (r) and dimension 

(r) for the position time series of a single oscillator of 
the FPU lattice for a range of different spatial resolutions 
r; see Ref. Q for further details. As can be appreciated 
in Fig. Ufa), it is not possible to obtain a m-independent 
scaling range for the correlation integral. Furthermore, 
from Fig. [2b) it is clear that the larger m the smaller 
is the scaling range for the correlation dimension, with 
remarkable statistical fluctuations at lower resolutions. 
Recalling that on large length scales the data cannot be 
distinguished from random noise [ill |23| ] , it becomes ev- 
ident that not even the dynamical character of the pro- 
vided signal can be inferred from this methodology. The 
problem can be traced back to the insufficient length of 
the employed time series. However, and since already 
for time series of 10 6 data points the computation time 
becomes prohibitively large, it is clear that the computa- 
tion of an adequate embedding dimension is impractical 
for a N = 32 FPU lattice. Thus, it is unfeasible to infer 
the dimensionality of this system from the available data 
by means of phase reconstruction methods 

The foregoing results impose severe restrictions to the 
computation of the LLE from time series, since it is 
known that the number of data points required to es- 
timate the LLE is about the square of that needed to 
estimate the embedding dimension (24|. Nevertheless we 
explore the extent of the aforementioned limitatio n by 
applying the algorithm proposed by Kantz in Ref. [25| . 
which does not depend explicitly on the knowledge of the 
correct embedding dimension, to our system. The results 
are displayed in Figs. [3]Ja) and (b), where S m (t), which 
is the temporal average of a suitable measure to the dis- 
tances between all neighboring trajectories to a reference 
trajectory inside a r-neighborhood for a given embedding 
dimension to, is plotted for low and high e values. It is 
clear that no temporal interval can be clearly identified, 
for all resolutions r studied, where S m (t) exhibits a linear 
increase with identical slope for any embedding dimen- 
sion value larger than m = 2. Therefore, no estimation of 
the LLE can be obtained from the employed time series, 
even in the strongly chaotic regime. It is important to 
mention that the data employed in the aforementioned 
calculation are exactly the same ones supplied to the 0-1 
test to compute the corresponding K value displayed in 
Fig. Q] and that these results are virtually unaltered when 
time series ten times larger are considered. Therefore, at 
least for this particular problem, the 0-1 test dramati- 
cally outperforms traditional phase space reconstruction 




FIG. 2: (Color online) (a) Correlation integral C m (r) and (b) 



dimension D\ ' (r) vs resolution r computed from the position 
time series of a single oscillator of a FPU lattice with N — 32 
and e = f 0, with embedding dimension m increasing from left 
to right in (a) and from bottom to top in (b). 




FIG. 3: (Color online) (a) Average local distance S m (t) vs 
time for the data set used in Fig. [2] with embedding dimen- 
sion increasing from bottom to top. Slope of dashed line cor- 
responds to the LLE value displayed in Fig. [Ha), (b) same 
as (a), but for a time series corresponding to e = 0.01. 



methods, which seemingly cannot cope with the dimen- 
sionality of the underlying dynamical system from where 
the time series originated. Our results are more conclu- 
sive than those in Ref . [4[ , where the LLE of the dissipa- 
tive eight-dimensional Lorenz 96 system could actually 
be computed with the method of Rosenstein et al. [26| 
and compared to the corresponding K values rendered 
by the 0-1 test. 



V. DISCUSSION AND CONCLUSIONS 

The preceding results clearly indicate that a quantita- 
tive estimate of the e T value corresponding to the SST 
can be efficiently obtained by means of the 0-1 test for 
chaos. The importance of this estimation stems from the 
fact that the SST has been detected in one-dimensional 
lattices with a (f> A interaction potential [lljj, as well as 
with Toda, smoothed Coulomb, and Lennard- Jones po- 
tentials [23| ; in an isotropic Heisenberg spin chain [28[ ; in 
a mean field XY chain [29j and in a coupled rotator chain 
which displays two thresholds separating two regions of 
weak chaos (occurring at low and high energies) from an 
intermediate region of strong chaos [3fJ] . It has also been 
detected in two- and three-dimensional lattices with two- 
well (f> 4 (3J HH and XY Heisenberg interactions [H, [34j . 
In particular, for a three-dimensional system consisting 



of a small number of ions confined in a Penning optical 
trap and forming a so-called microplasma the dependence 
of the LLE on the energy is very similar to that displayed 
by the aforementioned models where the detection of the 
SST has been performed [35j ] . Furthermore, a very de- 
tailed analysis of the behavior of the LLE as a function, 
not only of the energy but also of other parameters that 
control the geometry of the trap, is available [35j. 

At this point it is important to emphasize that the 
detection of the SST in the FPU model (or of a similar 
phenomenology in microplasmas, as already mentioned) 
does not require the precise knowledge of the LLE value: 
any observable (such as the asymptotic growth rate K) 
that shares the same dependence on the energy density 
with the LLE will be equally useful as well. From this ob- 
servation, and recalling that the trajectory of each indi- 
vidual ion in an optical trap can be observed and tracked 
experimentally, it is then highly feasible that much of the 
already known results obtained by the LLE for confined 
ions could be reproduced with less effort by means of 
the 0-1 test and corroborated by analogous calculations 
performed on experimental data. Furthermore , the max- 
imum ion number N = 40 employed in Ref. [35j | is of the 
same order than the oscillator number N = 32 in the 
present study; therefore, it can be reasonably expected 
that even the number of data points needed to charac- 
terize the microscopic dynamics of this microplasma by 
means of the 0-1 test for chaos could be of the same or- 
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der than that used to obtain the asymptotic growth rate 
K in Fig. [1] for the FPU lattice. Since, as already men- 
tioned, it is a difficult task to compute the LLE from 
experimental records, this simple and efficient technique 
could complement, and even become a viable alterna- 
tive when a rapid qualitative analysis is called for, to the 
approximate methods to the tangent space dynamics re- 
cently applied to characterize the microscopic dynamics 
of confined ions [36j . 

We wish to end our discussion by mentioning that 
"cuspy" patterns of the energy dependence of the LLE, 
in contrast to the mild transition displayed by the FPU 
model, show up in the presence of a thermodynamic 
phase transition, as in the case of the mean-field XY 
model and the two-dimensional lattice with (p 4 interac- 
tions previously mentioned. Now, it has been proposed 
that the appearance of singularities in the thermody- 
namic observables could be the effect of a suitable topo- 
logical transition in subspaces of configuration space [37} ■ 



It could be interesting to explore the possibility that in- 
formation about such detailed dynamical features could 
be probed by means of the herein presented methodol- 
ogy. In this context, another interesting possibility would 
be to test the feasibility of obtaining the abrupt tran- 
sitions just mentioned from time series stemming from 
variables different from the thermodynamic observables, 
which present singularities at phase transitions and are 
thus difficult to measure experimentally. 
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